function [tail] = tail_fit(varargin)
%
% tail current analsis. Double exponential fit to tail currents at a
% contant voltage, and boltzmanns fit to the amplitudes (and contant time)
% of the tails.
% Takes average double exp fit time constant and refits to get amplitudes.
% attempt to fit all at once leads to poor constraint and fitting (e.g.,
% hold tau's contant and fit all traces simultaneously).
%
% 11 July 2006 P. Manis
% Call from VCIV....
%
% Returns the tail structure, which contains:
% tail.fast = the parameters for the fast fit
% tail.slow = the parameters for teh slow fit
% tail.instant = the parameters for a direct measure (no fit).
%
global CURRENT VOLTAGE
global DFILE

tail = [];
time = make_time(DFILE);
RATES = DFILE.rate.*DFILE.nr_channel/1000;
hf = findobj('tag', 'tailfit');
if(isempty(hf))
    hf = figure('tag', 'tailfit');
else
    figure(hf);
    clf;
end;
figure(hf);
if(nargin == 0)
    tailstart = 210;
    taildur = 100;
    taildead = 3;
else
    tailstart = varargin{1};
    taildur = varargin{2};
    taildead = varargin{3};
end;
[a, ixt] =  min(abs(time(1,:)-tailstart)); % zero time for tail window measurement
[a, ixv0] = min(abs(time(1,:)-(tailstart-20))); % for voltage measurement
[a, ixv1] = min(abs(time(1,:)-(tailstart-10)));
[a, iinst] = min(abs(time(1,:)-(tailstart+5))); % time for "instantaneous" measurement
[a, ix1] =  min(abs(time(1,:)-(tailstart+taildur))); % duration of measurement
if(ix1 > size(CURRENT, 2))
    ix1 = size(CURRENT, 2);
end;
t0 = time(1,ixt);
[a, ix0] =  min(abs(time(1,:)-(tailstart+taildead))); % start time with deadwindow
tb = time(:,ix0:ix1)-t0;
vtail = mean(VOLTAGE(1,ixt:ix1));

ntr = size(CURRENT, 1);
npar = 3*ntr+2; % number of parameters in fit: tau 1, tau2, and then a0, a1 and offset for each
init = zeros(npar,1);
pmask = ones(npar, 1);
lbound = zeros(npar, 1);
ubound = zeros(npar, 1);
par1 = 1;
par2 = 1;
maxiter = 200;
% Smooth the current out a bit
for i=1:ntr
    fsamp = 1000/RATES(i); % get sampling frequency
    fco = 1000;		% cutoff frequency in Hz
    wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
    if(wco < 1) % if wco is > 1 then this is not a filter!
        [b, a] = fir1(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
        current_s(i,:) = filtfilt(b, a, CURRENT(i,:)); % filter all the traces...
    end
end
vm=zeros(ixv1-ixv0+1,1);
for j = 1:ntr
    % extract the voltage step information
    vm(j) = mean(VOLTAGE(j,ixv0:ixv1));
end;

avgflag = 0;
% to get initial estimate of time constants, we average the traces and fit
% one dual exp...
%
if(avgflag)
    avg = zeros(ix1-ix0+1,1);
    k = 0;
    for j = 1:ntr
        % extract the voltage step information
        vm(j) = mean(VOLTAGE(j,ixv0:ixv1));
        % now get an estimate of the time constants by fitting an averaged
        % trace.
        dz = 70;
        sa = current_s(j, ix0) - current_s(j, ix1);
        if(abs(sa) > dz)
            k = k + 1;
            if(k == 1)
                if(sa < 0)
                    avg = -current_s(j,ix0:ix1);
                else
                    avg = current_s(j,ix0:ix1);
                end;
            elseif(sa < 0)
                avg = avg - current_s(j,ix0:ix1);
            else
                avg = avg + current_s(j,ix0:ix1);
            end;
        end;

    end;
    avg = avg/k;
    [FP, CHISQ, NITER] = mrqfit('exponential', [avg(end), avg(1)-avg(end), 5, 25, 50], tb(1,:), avg, ...
        [], [], [-500, -1000, 0.5, -1000, 10], [500, 1000, 20, 1000, 80], 200, 1e-3 );
    t1a = FP(3);
    t2a = FP(5);
    plot(tb(1,:), avg, '-k');
    yt = FP(1)+FP(2)*exp(-tb(1,:)/FP(3))+FP(4)*exp(-tb(1,:)/FP(5));
    hold on;
%    plot(tb(1,:), yt, '-r');
    pause(3);
    clf;
else
    t1a=5;
    t2a=40;
end;
%ha=subplot('Position', [0.1 0.5 0.8 0.4]);
%plot(ha, tb', current_s(:,ix0:ix1)', '-k');
%hold on;

clf;
ha=subplot('Position', [0.1 0.5 0.8 0.4]);
t1=[];
t2=[];
for i = 1:ntr
    a0 = current_s(i, ix1);
    a1 = current_s(i, ix0)-a0;
    [FP, CHISQ, NITER] = mrqfit('exponential', ...
        [a0, a1*0.7, t1a, a1*0.3, t2a], tb(i,:), current_s(i,ix0:ix1), ...
        [], [], [-500, -1000, 0.5, -1000, 10], [500, 1000, 20, 1000, 100], 200, 1e-3 );
    t1(i) = FP(3);
    t2(i) = FP(5);
    j = i*3;
    init(j) = FP(1);
    init(j+1) = FP(2);
    init(j+2) = FP(4);
    plot(ha, tb(1,:), current_s(i,ix0:ix1), '-k');
    yt = FP(1)+FP(2)*exp(-tb(1,:)/FP(3))+FP(4)*exp(-tb(1,:)/FP(5));
    hold on;
    plot(ha, tb(1,:), yt, '-r');
end;
%pause(2);

init(1) = mean(t1);
init(2) = mean(t2);

if1 = [1:ntr];
maxiter = 500;
tryallfit = 0;
if(tryallfit)
    lbound(1) = init(1)-std(t1)*2; ubound(1) = init(1)+std(t1)*2;
    lbound(2) = init(2)-std(t2)*2; ubound(2) = init(2)+std(t2)*2;
    %pmask(1) = 0; pmask(2) = 0;
    % ar1 = FP(2)/(FP(4)+FP(2)+FP(1));
    % ar2 = FP(4)/(FP(4)+FP(2)+FP(1));
    for i = 1:ntr
        j = i*3;
        %     init(j) = current_s(i,ix1);
        %     init(j+1) = ar1*(current_s(i,ix0)-init(j));
        %     init(j+2) = ar2*(current_s(i,ix0)-init(j));
        lbound(j) = init(j)*0.5; lbound(j+1) = init(j+1)*0.5; lbound(j+2) = init(j+2)*0.5;
        ubound(j) = init(j)*2; ubound(j+1) = init(j+1)*2; ubound(j+2) = init(j+2)*2;
    end;
    % initial fit is done on 4 traces. These parameters are then used
    % to seed the final fit.
    %function [fit_error,pars]=curve_fitting(x_data,y_data,algorithm,search,model,...
    %   order,init,PMASK,LBOUND,UBOUND,Par_1,Par_2, maxiter)

    %[1,floor(0.27*ntr),floor(0.45*ntr),floor(0.66*ntr),ntr];
    [f, pars] = curve_fitting(tb(if1,:), current_s(if1,ix0:ix1), 'levenberg', 'cubic', 74, ...
        npar, init, pmask, lbound, ubound, par1, par2, maxiter)


    a1=[];
    a2=[];
    y=[];
    i1 = 0;
    for i = if1
        i1 = i1 + 1;
        j = i * 3;
        y(i1,:) = pars(j)+(pars(j+1).*exp(-tb(i,:)/pars(1)))+(pars(j+2).*exp(-tb(i,:)/pars(2)));
        a1(i1) = pars(j+1);
        a2(i1) = pars(j+2);
    end;

    plot(ha, tb', y', '-g');
else
    t1a=init(1);
    t2a=init(2);
    a1=[]; a2=[]; a0=[];
    for i = 1:ntr
        a0i = current_s(i, ix1);
        a1i = current_s(i, ix0)-a0i;
        [FP, CHISQ, NITER] = mrqfit('exponential', ...
            [a0i, a1i*0.7, t1a, a1i*0.3, t2a], tb(i,:), current_s(i,ix0:ix1), ...
            [], [1 1 0 1 0], [-500, -1000, 0.5, -1000, 10], [500, 1000, 20, 1000, 100], 200, 1e-3 );
        t1(i) = FP(3);
        t2(i) = FP(5);
        j = i*3;
        init(j) = FP(1);
        init(j+1) = FP(2);
        init(j+2) = FP(4);
        a1(i) = FP(2);
        a2(i) = FP(4);
        a0(i) = FP(1);
        %  plot(ha, tb(1,:), current_s(i,ix0:ix1), '-k');
        yt = FP(1)+FP(2)*exp(-tb(1,:)/FP(3))+FP(4)*exp(-tb(1,:)/FP(5));
        hold on;
        plot(ha, tb(1,:), yt, '-g');
    end;
    pause(5);
end;

hb=subplot('Position', [0.1 0.1 0.35 0.35]);

%[A0 A1 K1 X1];
% MRQFIT(F, P, X, Y, SIG, VP, LB, UB, IMAX, TOL)
[FP2, CHISQ2, NITER2] = mrqfit('boltzmann', [min(a2), max(a2)-min(a2), 1, -45], vm(if1), a2, ...
    [], [], [-1000, -1000, 0.01, -120], [1000, 1000, 25, +50], 200, 1e-5 );
yb2 = FP2(1) + (FP2(2)-FP2(1))./(1+exp(-(vm-FP2(4))*FP2(3)));

[FP1, CHISQ1, NITER1] = mrqfit('boltzmann', [min(a1), max(a1)-min(a1), 1, -45], vm(if1), a1, ...
    [], [], [-1000, -1000, 0.01, -120], [1000, 1000, 25, +50], 200, 1e-5 );
yb1 = FP1(1) + (FP1(2)-FP1(1))./(1+exp(-(vm-FP1(4))*FP1(3)));

[FP3, CHISQ3, NITER3] = mrqfit('boltzmann', [min(a2), max(a2)-min(a2), 1, -45], vm(if1), current_s(:,iinst), ...
    [], [], [-1000, -1000, 0.01, -120], [1000, 1000, 25, +50], 200, 1e-5 );
yb3 = FP3(1) + (FP3(2)-FP3(1))./(1+exp(-(vm-FP3(4))*FP3(3)));
plot(hb, vm(if1), a1, 'go');
hold on;
plot(hb, vm, yb1, 'g-');
plot(hb, vm(if1), a2, 'rx');
plot(hb, vm, yb2, 'r-');
plot(hb, vm(if1), current_s(if1,iinst), 'bs');
plot(hb, vm, yb3, 'b-');
legend('fast', 'fast-fit', 'slow', 'slowfit', 'instantaneous', 'instfit');
fprintf(1, ['t1 = %7.4f Vh1 = %7.2f  k1 = %7.4f A1 = %7.2f (green)\n', ...
    't2 = %7.4f Vh2 = %7.2f  k2 = %7.4f A2 = %7.2f (red)\n', ...
    'inst = %7.4f Vh3 = %7.2f  k3 = %7.4f A3 = %7.2f (blue)\n']',...
    t1a, FP1(4), 1/FP1(3), FP1(2)-FP1(1), ...
    t2a, FP2(4), 1/FP2(3), FP2(2)-FP2(1), ...
    time(1,iinst), FP3(4), 1/FP3(3), FP3(2)-FP3(1));

tail.bfast = FP1;
tail.bslow = FP2;
tail.binstant = FP3;
tail.taufast = t1a;
tail.tauslow = t2a;
tail.vm=vm;
tail.a1 = a1;
tail.a2 = a2;
tail.inst = current_s(:,iinst);
tail.vtail = vtail;
